A Family of Equations of State Based on Lattice QCD: 
Impact on Flow in Ultrarelativistic Heavy-Ion Collisions 
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We construct a family of equations of state within a quasiparticle model by relating pressure, 
energy density, baryon density and susceptibilities adjusted to first-principles lattice QCD calcula- 
tions. The relation between pressure and energy density from lattice QCD is surprisingly insensitive 
to details of the simulations. Effects from different lattice actions, quark masses and lattice spacings 
used in the simulations show up mostly in the quark-hadron phase transition region which we bridge 
over by a set of interpolations to a hadron resonance gas equation of state. Within our optimized 
quasiparticle model we then examine the equation of state along isentropic expansion trajectories at 
small net baryon densities, as relevant for experiments and hydrodynamic simulations at RHIC and 
LHC energies. We illustrate its impact on azimuthal flow anisotropies and transverse momentum 
spectra of various hadron species. 
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I. INTRODUCTION 

In the last few years, much evidence has been accumu- 
lated for the applicability of hydrodynamics in describing 
the expansion stage of strongly interacting matter cre- 
ated in relativistic heavy- ion collisions [1~7]. Hydrody- 
namics describes the collective flow of bulk matter from 
an initial state just after reaching thermalization up to 
the kinetic freeze-out stage. The heart of hydrodynamics 
is the equation of state (EoS) which relates thermody- 
namically the pressure p of the medium to its energy 
density e and net baryon density ns (or, equivalently, 
to its temperature T and baryon chemical potential /is). 
Specifically, the parameter controlling the acceleration of 
the fluid, i.e. the build-up of collective flow, by pressure 
gradients in the system is the speed of sound, given by 

„2 _ dp 
— de- 

While most existing hydrodynamic simulations have 
used a realistic hadron resonance gas EoS below the de- 
conflnement transition (either with full [1, 4, 5] or par- 
tial [2, 8-11] chemical equilibrium among the hadron 
species), they have usually relied on simple analyt- 
ical models for the EoS of the quark-gluon plasma 
(QGP) above the transition, based on the assumption 
of weak coupling among the deconflned quarks and glu- 
ons. This assumption is, however, inconsistent with the 
phenomenological success of hydrodynamics which re- 
quires rapid thermalization of the QGP [12] and therefore 
strong interactions among its constituents [13-16]. In- 
deed, lattice QCD calculations of the QGP pressure and 
energy density show that they deviate from the Stefan- 
Boltzmann limit for an ideal gas of non-interacting 
quarks and giuons even at temperatures T > 3 Tc (with 
Tc as pseudo-critical temperature), by about 15-20% [17- 
19]. Miraculously, however, the deviations are of similar 
magnitude in both p and e such that, for T > 2Tc, the 



squared speed of sound = ^ ~ ^ [19], just as expected 
for a non-interacting gas of massless partons. In spite of 
the evidence for strong interactions among the quarks 
and giuons in the QGP seen in both p{T) and e(T), the 
stiffness and accelerating power of the lattice QCD equa- 
tion of state is thus indististinguishable from that of an 
ideal parton gas (at least for temperatures T > 2Tc), 
such as the one used above Tc in most hydrodynamical 
simulations. 

On the other hand, at T < 2 Tc the speed of sound 
extracted from lattice QCD drops below the ideal gas 
value Cs = l/"\/3, reaching a value that is about a fac- 
tor of 3 smaller near Tc [19]. This leads to a signiflcant 
softening of the QGP equation of state relative to that of 
an ideal massless gas exactly in the temperature region 
Tc<T <2Tc explored during the early stages of Au-|-Au 
collisions at RHIC [1, 2, 4, 5, 8]. To explore the sensi- 
tivity of the flow pattern seen in the RHIC data to such 
details of the EoS near the quark-hadron phase transi- 
tion, the hydrodynamic evolution codes must be supplied 
with an EoS that faithfully reproduces the lattice QCD 
results above Tc- To construct such an EoS, and to test 
its influence on the collective flow generated in RHIC and 
LHC collisions, are the main goals of this paper. 

Our approach is based on the quasiparticle model [20- 
29] which expresses the thermodynamic quantities as 
standard phase space integrals over thermal distribu- 
tion functions for quasiparticles with medium dependent 
properties. In the present paper we follow the philoso- 
phy [20-28] that the interaction effects in the QGP can be 
absorbed into the quasiparticle masses and a vacuum en- 
ergy all of which depend on the temperature and baryon 
chemical potential. This is known to produce good fits 
to the lattice QCD data both at vanishing [20-23] and 
non-vanishing [25-27] baryon chemical potential. How- 
ever, since this approach uses on-shell spectral functions 
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for tfie quasiparticlcs, it implicitly assumes zero resid- 
ual interactions (i.e. infinite mean free paths) for them, 
which is inconsistent with the low viscosity and almost 
ideal fluid dynamical behaviour of the QGP observed at 
RHIC. Peshier and Gassing [29] have shown that it is pos- 
sible to generalize the quasiparticle description to include 
a finite (even large) coUisional width in the spectral func- 
tions, without significantly affecting the cjuality of the 
model fit to the lattice QCD data for the EoS at /i^ = 0. 
Since hydrodynamics only cares about the EoS, but not 
about its microscopic interpretation, we here opt for the 
simpler, but equally successful approach using on-shell 
quasiparticlcs to fit the lattice QCD EoS. 



The quasiparticle EoS for the QGP above Tc does 
not automatically match smoothly with the hadron reso- 
nance gas EoS below Tc- Although the gap between the 
two branches of the EoS is much smaller here than for 
the previously used models which assume non-interacting 
quarks and gluons above Tc [1-5, 8-11], a certain de- 
gree of ambiguity remains in the interpolation process. 
We explore a set of different interpolation presc;riptions, 
yielding a family of equations of state which exhibit slight 
differences in the phase transition region, and study their 
dynamical consequences. 



Our paper is organized as follows: In Sec. II we show 
that our quasiparticle model provides an efficient and ac- 
curate parametrization of lattice QCD results for Nf = 2 
flavors both at /is = and /i^ 7^ 0. We also extract 
the isentropic expansion trajectories followed by fully 
equilibrated systems. In that Section, the quasiparticle 
parametrization is continued below To, down to temper- 
atures of about 0.75 Tc where the lattice QCD data end. 
In Sec. Ill we proceed to the physically relevant case of 
Nf =2 + 1 flavors and furthermore match the quasiparti- 
cle EoS above Tc to a hadron resonance gas EoS below Tc. 
Variations in the matching procedure lead to a family of 
equations of state with slightly different properties near 
Tc- The transition to a realistic hadron resonance gas 
picture below Tc means that these EoS can now be used 
down to much lower temperatures to make explicit con- 
tact with the experimentally observed final state hadrons 
after decoupling from the expanding fluid. In Sec. IV we 
use this family of EoS for hydrodynamic calculations of 
the differential elliptic flow V2{pt) for several hadronic 
species in Au-|-Au collisions at the top RHIC energy and 
compare with experimental data. We find some sensitiv- 
ity to the details of the interpolation scheme near Tc, as 
long as an EoS is used that agrees with the lattice QCD 
data for energy densities e > 4GeV/fm^. We conclude 
that Section with a few predictions for Pb+Pb collisions 
at the LHC. A short summary is presented in Sec. V. 



II. QUASIPARTICLE DESCRIPTION OF THE 
EOS FROM LATTICE QCD FOR Nf = 2 

A. The quasipeirticle model 

Over the years, several versions of quasiparticle models 

have been developed to describe lattice QCD data for 
the QCD equation of state [20-27, 29]. They differ in 
the choice and number of parameters and in the details 
of the underlying microscopic picture but generally yield 
fits to the lattice QCD data which are of similar quality. 
In this subsection we quickly review the essentials of the 
model described in [20] which will be used here. 

In our quasiparticle approach the thermodynamic pres- 
sure is written as a sum of contributions associated with 
medium modified light quarks q, strange quarks s, and 
gluons g [20]: 

p{T,{fla})= Y^Pa-B{T,{fia}), (1) 
a=q,s,g 

with partial pressures 

Here = (exp([tJa T Ma]/T) + Sa)^^ are thermal 
equilibrium distributions for particles and antiparticles, 
with Sq^s = 1 for fermions and Sg = —1 for bosons. 
da represents the spin-color degeneracy factors, with 
dq = 2NqNc = 12 for the Ng = 2 light quasi-quarks, 
ds = 2Nc = 6 for the strange quasi-quarks, and dg = 
— 1 = 8 for the right-handed transversal quasi-gluons 
(with the left-handed ones counted as their antiparticles). 
Since the pressure integral in Eq. (2) is dominated by 
thermal momenta of order k ^ T, weak coupling pertur- 
bation theory suggests [30, 31] that the dominant propa- 
gating modes are transversal plasmons with gluon quan- 
tum numbers {g) and quark-like excitations, whereas lon- 
gitudinal plasmons are exponentially suppressed. Our 
model assumes that this remains true near Tc where per- 
turbation theory is not expected to be valid. 

We are interested in the application of this EoS to 
heavy ion collisions where strangeness is conserved at its 
initial zero value, due to the very short available time. 
This strangeness neutrality constraint allows to set fXg = 

0. The isospin chemical potential /i/ = (/i^ — ^-id)/'2 is 
fixed by the net electric charge density of the medium; we 
assume zero net charge of the fireball matter created near 
midrapidity at RHIC as well as equal masses for the u and 
d quasi-quarks such that /i/ = and we have only one 
independent chemical potential ^„ = /id = = /is/S 
where /^s is the baryon number chemical potential. 

The quasiparticles are assumed to propagate on-shell, 

1. e, with real energies uja given by dispersion relations 
of the type u)a = ^/W+Tn^J^\Jlq) , known to hold for 
weakly interacting quarks and gluons with thermal mo- 
menta k ^ T. Again the model assumes that this struc- 
ture holds true also near Tc where perturbation theory 
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presumably breaks down. In order to directly compare 
our Quasiparticle Model (QPM) with lattice QCD re- 
sults, we include nonzero bare quark masses Wao and 
adjust them to the values used in the lattice simula- 
tions through ml = m^g + Ua [32] where Ha denotes 
the self energy. For gluonic modes we use nigo = 0. 
For Ua we employ an ansatz inspired by the asymptotic 
form of the gauge independent hard thermal/dense loop 
(HTL/HDL) self-energies which depend on T, Hg, rUao, 
and the running coupling g"^ as follows [30, 32] : 



n„ 



n„ 




(3) 



(4) 



(5) 



The jUa in the dispersion relations thus denote effec- 
tive quasiparticle masses due to the dynamically gener- 
ated self-energies Ua- The mean field interaction term 
B{T,iiq) in Eq. (1) is determined by thermodynamic 
self-consistency and stationarity of the thermodynamic 
potential under functional variation of the self-energies, 
Sp/SHa = [33]. As a consequence, B{T,^q) is evalu- 
ated in terms of an appropriate line integral in the T-^q 
plane, with integration constant B(Tc) adjusted to the 
lattice results [20]. 

All other thermodynamic quantities follow straight- 
forwardly from the stationarity condition and standard 
thermodynamic relations. For example, the entropy den- 
sity reads 



a=q,s,g 



(6) 



So, = 



da 

27r2 



Jo 



k^dk 



^ rp a ~^fa ) rp if a fa^ 



n„ 



while the net quark number density Uq = Sub is given 
through 

-/,-)• (7) 

Although the form of our ansatz for the quasiparticle 
masses (i.e. the specific interplay between the param- 
eters mao,T, and /Xg) is inspired by perturbation the- 
ory, our model becomes non-perturbative by replacing 
the perturbative expression for the running coupling 
in Eqs. (3-5) by an effective coupling whose depen- 
dence on T and jiq is parametrized and fitted to the non- 
pcrturbativc (T, /iq)-dcpendcncc of the thermodynamic 
functions from lattice QCD. The (T, /Xg)-dependence of 
is constrained by Maxwell's relation for p which takes 
the form of a quasi-linear partial differential equation 



dllq 



dG^ , 



(8) 



here o^^ , ay and b depend on T, /ig and G^ (sec 
Refs. [20, 34] for details). This flow equation is solved by 
the method of characteristics, starting from initial con- 
ditions on a Cauchy surface in the T-iJ.q plane. One pos- 
sibility is to parameterize G^ at /U, = such that lattice 
QCD results for vanishing quark chemical potential are 
reproduced, and to use the flow equation for extrapola- 
tion to non-zero yUg. As a convenient parametrization of 
G^{T,iiq=0) we find [35] 



G^{T,f,q=0) 



'-'2~loop(^c. 



(9) 



Here, in order to recover perturbation theory in the high 
temperature limit, G2_ioop is taken to have the same form 
as the perturbative running coupling at 2-loop order: 



<-^2-loop(^) 



1 



2/3i log(loge2) 



(10) 



with /3o = i(lliVc - 2Nj) and /3i = ^{SiN^ - ISNjNc + 
SNf/N^). The scale ^ is parametrized phenomenologi- 
cally as ^ = A(T — Ts)/Tc, with a scale parameter A and 
a temperature shift Tg which regulates the infrared di- 
vergence of the running coupling by shifting it somewhat 
below the critical temperature Tc- Below the phase tran- 
sition, we postulate a continuous linear behavior of the 
effective coupling. The parametrization (9,10) turns out 
to be flexible enough to describe the lattice QCD results 
accurately down to about T « 0.75 Tg. In contrast, using 
a pure 1-loop or 2-loop perturbative coupling together 
with a more complete description of the plasmon term 
and Landau damping restricts the quasiparticle approach 
to T > 2Tc [36]. [Similar quality fits can be achieved 
in that approach, without giving up its more accurate 
form of the HTL/HDL self energies, by adopting a simi- 
lar non-perturbative modification of the running coupling 
as adopted here [23].] 

The model described in this subsection was successfully 
applied to QCD lattice data in the pure gauge sector in 
Ref. [20], and to first lattice QCD calculations at /x, ^ 
in Ref. [37]. In the following subsection we test it on 
recent lattice QCD data for Nf = 2 dynamical quark 
flavors at zero and non-zero /iq, and in the next section 
we consider the realistic case of Nf = 2 -|- 1 flavors with 
the aim of providing an EoS suitable for hydrodynamic 
simulations of heavy-ion collisions. 



B. Thermodynamics of Nf — 2 quark flavors 

We begin with the case of Nf = 2 dynamical quark 
flavors at zero quark chemical potential and confront the 
QPM with lattice QCD results obtained by the Bielefeld- 
Swansea collaboration [17]. These simulations were per- 
formed with temperature dependent bare quark masses 
mao{T) = XaT where Xg = Q and Xq = 0.4 [17]. For 
Nf = 2 light quark flavors we can set = in the 
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QPM expressions. Fig. 1 shows the lattice QCD data 
for the scaled pressure p{T)/T'^ together with the QPM 
fit; the fit parameters given in the caption were obtained 
by the procedure described in Ref. [38]. The raw lattice 
data were extrapolated to the continuum by multiplying 
the pressure in the region T > Tc by a constant factor 
d = 1.1, following an estimate given in [17, 40] who ad- 
vocate a range of 10-20% due to finite size and cutoff 
effects. (Note that this estimated correction factor does 
not necessarily have to be independent of T, as assumed 
here.) 




T/T, 



FIG. 1: (Color online) Comparison of the QPM with lat- 
tice QCD results (symbols) for the scaled pressure p/T* as 
a function of T/Tc for Nf = 2 and fiq = 0. The raw lat- 
tice QCD data from [17] have been continuum extrapolated 
as described in the text. The QPM parameters are A — 4.4, 
Ts = 0.67n, b = 344.4 and B{Tc) = O.SIT*, with Tc = 175 
MeV as suggested in [39]. The horizontal line indicates the 
Stefan-Boltzmann value psb/T" = cq = (32 -I- 21iV/)7r^/180 
for Nf = 2. 

Having demonstrated the ability of the QPM to suc- 
cessfully reproduce lattice EoS data along the — 
axis, we can now exploit recent progress in lattice QCD 
with small non-vanishing chemical potential to test its 
ability to correctly predict the thermodynamic functions 
at non-zero fj,q. In Ref. [41] finite- /i^ effects were eval- 
uated by expanding the pressure into a Taylor series in 
powers of (fiq/T) around fiq — 0, 



n=0,2,4,. 



Cn{T) 



(11) 



where co(T) = p{T,Hq=0)/T^ is the scaled pressure 
at vanishing quark chemical potential. The coefficients 
C2(r), C4(r), cg{T) were extracted from the lattice by 
numerically evaluating appropriate /K^-derivatives of the 
logarithm of the partition function InZ = pV/T [41], viz. 



c„(T) 



1 d''ip/T* 



n! difiq/ry 



(12) 



Note that computing the coefficients c„, n > 2, from 
these expressions is easier on the lattice than determining 
the pressure at /is = 0, co{T), since the latter requires 
an integration over T and a separate lattice simulation 
at T = 0. For this reason Ref. [41] has no results for 
co(T). Since the simulations in Ref. [41] were done with 
different parameters than those analyzed in Fig. 1 [17], it 
is not immediately clear that the QPM parameters fitted 
to the results of Ref. [17] can also be used to describe the 
simulations reported in [41]. When analyzing the lattice 
data of [41] we therefore refit the QPM parameters to 
the lattice results for C2(T) (see dashed line and squares 
in Fig. 4 below) and then assess the quality of the model 
fit by its ability to also reproduce C4{T) and Ce{T) ex- 
tracted from the same set of simulations, as well as other 
thermodynamic quantities calculated from these coeffi- 
cients through Taylor expansions of the type (11). The 
QPM parameters obtained by fitting C2(T) from [41] are 
[35] A = 12.0, T, = 0.87 Tc, and b = 426.05 (again using 
Tc = 175 MeV) [42]. 

Evaluation of the derivatives in (12) within the QPM 
is straightforward; for explicit analytical expressions for 
•^2,4,6(7") we refer the reader to equations (6, 7,8) in the 
second paper of Ref. [35] . That paper also shows that the 
quasiparticle model gives an excellent fit to C2(T) from 
[41], and that with the same set of parameters the QPM 
expressions for c^iT) and ce{T) yield impressive agree- 
ment with the lattice data [41], too. In particular, several 
pronounced structures seen in C4(T) and cg{T) are quan- 
titatively reproduced [35]. This constitutes a stringent 
test of the efficiency of our QPM parametrization. 

We here use these first three expansion coefficients 
£2,4,6(2^) to write down truncated expansions for the net 
baryon density ns = dp/dfiB and the corresponding 
baryon number susceptibility xb — duB/d^B which is 
a measure of fluctuations in : 



nB{T,^B) ^ 

3 ^ ' \3TJ 3 

XB{T,fiB) 



(13) 



C4(T)(g)V2ca(T)(^V 



.3r/ 



y2 



(14) 



+ (I?)' + y-(-'(f?)' 



These yield a truncated result for p{T,fiq). 



In Fig. 2, the truncated QPM results for ub/T^ and 
Xb/T"^ are compared for various values of ^b/Tc with 
lattice QCD results that were obtained from Eqs. (13) 
and (14) with the coefficients C2,4,6(2^) from [41]. We find 
good agreement with the lattice results; even below 
where our QPM parametrization is not well justified and 
should be replaced by a realistic hadron resonance gas 
(see Sec. Ill), the deviations are small but increase with 
increasing ^b/Tc- All in all, the QPM model appears to 
provide an efficient and economic parametrization of the 
lattice data down to T '--^ 0.75 Tc. 

Within the QPM model we can assess the truncation 
error made in Eqs. (13) by comparing this expression 
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T/T, 

FIG. 2: (Color online) Scaled baryon density ns/T'^ (upper 
panel) and baryon number susceptibility xb/T"^ (lower panel) 
as a function of T/Tc, for fis/T^^ 2.4, 1.8, 1.2, 0.6 (from top 
to bottom). QPM results from the truncated expansions (13) 
and (14) (solid lines) are compared with lattice QCD data 
(symbols) from [41] for Nf — 2. Dashed lines in the upper 
panel represent the full QPM result (7) for ub = riq/S. The 
QPM parameters are A = 12.0, Ts = 0.87 Tc, and b = 426.05, 
for Tc = 175 MeV. 

with the exact result (7) (dashed lines in the upper panel 
of Fig. 2). The authors of [41] estimated the error in- 
duced in Eq. (11) by keeping only terms up to n = 4 to 
remain < 10% for fis/T < 3. Here we keep the terms 
~ {fis/T)^ and, as the upper panel of Fig. 2 shows, the 
resulting truncated expessions for the baryon density ns 
agree with the exact results within the linewidth as long 
as ^b/Tc < 1.8. For hb/Tc ~ 2.4 we see significant devi- 
ations between the truncated and exact expressions near 
T = Tc which, however, can be traced back to an arti- 
ficial mechanical instability dp/dns < induced by the 
truncation. Similar truncation effects near T = Tc are 
stronger and more visible in the susceptibility xb (lower 
panel of Fig. 2). In both cases the full QPM expression 
is free of this artifact and provides a thermodynamically 
consistent description. 

We next compare the Taylor series expansion coef- 
ficients of the energy and entropy densities given in 



Ref. [39] with our model. We have the following decom- 
positions [39]: 



oo 



(15) 



n=0 



s = s{T, fi.^O) + ^ ((4-n)c„(r) 



c;(T))(^ 



with p from (11), c[^{T) = Tdcn{T)/dT, and 
s{T,pi,=0) = T^4co{T) + c',{T)). 



(16) 



Since these expressions contain both c„(r) and their 
derivatives with respect to T, c^(T), they provide a more 
sensitive test of the model than considering the pressure 
alone. The expressions (15) can be read as Taylor series 
expansions with expansion coefficients 

^ = ^e„(T)(^)", e„(T) = 3c„(T) + c'„(T), (17) 

s„(T) = (4-n)c„(r) + cUT). 



s 



n 



Figure 3 shows a comparison of the QPM results for 62,4 
and S2,4 (obtained through fine but finite difference ap- 
proximations of the c„(T)) with the corresponding lattice 
QCD results from Ref. [39]. The QPM parameters are 
the same as in Fig. 2, and the agreement with the lattice 
data is fairly good. The pronounced structures observed 
in the vicinity of the transition temperature are a result 
of the change in curvature of G'^(T, /ig=0) at T = Tc (see 
Eq. (9)). Note that the derivatives c'^{T) were estimated 
in [39] by finite difference approximations of the available 
lattice QCD results for c„(T). After adjusting the differ- 
ence approximation in our QPM to the lattice procedure, 
tlie pronounced structures in the vicinity of Tc are much 
better reproduced [43]. 

We close this subsection with a calculation of the quark 
number susceptibilities which play a role in the calcula- 
tion of event-by-event fluctuations of conserved quanti- 
ties such as net baryon number, isospin or electric charge 
[44-47]. Across the quark- hadron phase transition they 
are expected to become large. For instance, the peak 
structure in C4(r) (which for small fJ-s/T gives the dom- 
inant /^B-dependence of xb, see Eq. (14)) can be in- 
terpreted as an indication for critical behavior. Quark 
number susceptibilities have been evaluated in lattice 
QCD simulations by Gavai and Gupta [48], using con- 
stant bare quark masses m^o — 0.1 Tc with Tc fixed by 
nip/Tc — 5.4. Introducing separate chemical potentials 
for u and d quarks and considering a simultaneous ex- 
pansion of the QCD partition function Z{T, ii-^, iJ-d) in 
terms of /x„ and /i^j, the leading /Zu^d-iiidependent contri- 
bution to the quark number susceptibility Xq = ^Xb can 
be expressed in terms of Xuu , Xud and Xdd where 



Xab 



(18) 
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FIG. 3: (Color online) Comparison of the Taylor series expan- 
sion coefficients for e„ {T) (squares/dashed black lines) and 
s„{T) (circles/solid red lines) for Nf = 2 from [39] with the 
QPM (same parameters as in Fig. 2). [Upper panel: n = 2. 
Lower panel: n = 4.] For details see text. 



These linear quark number susceptibilities can be related 
to the Taylor series expansions in (11) and (14) through 

C2{T) ^ iXuu + "^Xud + Xdd) ■ (19) 

For m„ = rrid one finds Xuu = Xdd- In Fig- 4 we compare 
lattice QCD results [48] for {xuu+Xud)/T'^ = C2{T) with 
a QPM fit. The QPM parameters are adjusted to the 
lattice data from [48], after extrapolating the latter to 
the continuum by multiplying with a factor d = 0.465 as 
advocated in [49]. For comparison, we also show C2{T) 
from [41] and the corresponding QPM parametrization 
from Fig. 2. Note that the latter data have not yet been 
extrapolated to the continuum. If we performed a contin- 
uum extrapolation of the C2(T) data from [41] by a factor 
d = 1.1 for T > Tc as in the case of Co(T) (cf. Fig. 1), 
both results would agree at large T within 1%. In the 
transition region some deviations would remain, due to 
the different bare quark masses and actions employed in 
Refs. [41] and [48]. 



1 




FIG. 4: (Color online) Comparison of the QPM result for 
iXuu + Xud)/T^ (solid line) with lattice QCD data (cir- 
cles) from [48] for Nf = 2, extrapolated to the continuum 
as suggested in [49]. The QPM parameters are A — 7.0, 
Ts = 0.76 Tc, and 6 = 431, with Tc = 175 MeV. For com- 
parison, we also show lattice QCD data for C2{T) for Nf — 2 
from [41] (squares) together with the the corresponding QPM 
fit (dashed line), using the same parameters as in Fig. 2. 



C. Isentropic trajectories for Nf = 2 quark flavors 

Ideal relativistic hydrodynamics [1~6] is considered to 
be the appropriate framework for describing the expan- 
sion of strongly interacting quark-gluon matter created in 
relativistic heavy-ion collisions. This approach requires 
approximate local thermal equilibrium and small dissipa- 
tive effects. Since the fireballs created in heavy-ion exper- 
iments are small, pressure gradients are big and expan- 
sion rates are large, thermalization must be maintained 
by sufRciently fast momentum transfer rates resulting in 
microscopic thermalization time scales which are short 
compared to the macroscopic expansion time. The hy- 
drodynamic description remains valid as long as the par- 
ticles' mean free paths are much smaller than both the 
geometric size of the expanding fireball and its Hubble 
radius. 

The hydrodynamic equations of motion result from 
the local conservation laws for energy-momentum and 
conserved charges, d^j_T^^''{x) — and 9^jf (x) — 0. Here, 
T'^" denotes the energy-momentum stress tensor and 
j,f the four-current of conserved charge i at space-time 
coordinate x. Heavy-ion collisions are controlled by 
the strong interaction which conserves baryon number, 
isospin, and strangeness. If we assume zero net isospin 
and strangeness densities in the initial state, only the 
conservation of the baryon number four-current needs 
to be taken into account dynamically. 

The ideal fluid equations are obtained by assuming 
locally thermalized momentum distributions in which 
case r^'' and take on the simple ideal fluid forms 
T'^'^ = {e + p)u''u'' -pg''" and = n^u'^ [50]. Here gt"" 
is the Minkowski metric, u^{x) the local four- velocity of 
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the fluid, and e{x), p{x) and nsix) denote the energy 
density, pressure, and baryon density in the local fluid 
rest frame. The resulting set of 5 equations of motion for 
6 unknown functions is closed by the EoS which relates 
p, e, and ub- This is where the lattice QCD data and 
our QPM parametrization of the lattice EoS enter the 
description of heavy-ion collision dynamics. 

Once the initial conditions are specified, the further dy- 
namical evolution of the collision fireball is entirely con- 
trolled by this EoS. Specifically, the accelerating power of 
the fluid (i.e. its reaction to pressure gradients which pro- 
vide the thermodynamic force driving the expansion) is 
entirely controlled by the (temperature dependent) speed 
of sound, Cg = \fdplde. To the extent that ideal fluid 
dynamics is a valid description and/or dissipative effects 
can be controlled, the observation of collective flow pat- 
terns in heavy-ion collisions can thus provide constraints 
on the EoS of the matter formed in these collisions. 

Ideal fluid dynamics is entropy conserving, i.e. the 
specific entropy a = s/ub of each fluid cell (where s is 
the local entropy density) stays constant in its comov- 
ing frame. Although different cells usually start out with 
different initial specific entropies, and thus the expand- 
ing fireball as a whole maps out a broad band of widely 
varying s/ns values, each fluid cell follows a single line of 
constant s/ub in the T—fiB phase diagram. It is there- 
fore of interest to study the characteristics of such isen- 
tropic expansion trajectories, in particular the behavior 
of p/e or = ^ along them. 

The isentropic trajectories for different values oi s/ub 
follow directly from the flrst principles evaluation of the 
lattice EoS and its QPM parametrization considered in 
the previous subsection. For Nf = 2 dynamical quark fla- 
vors, the truncated Taylor series expansions for baryon 
number and entropy density with expansion coefficients 
c„(T) and s„(r) according to (17) were employed in Ref. 
[39] to determine the isentropic trajectories for s/ub = 
300, 45, 30, sampling those regions of the phase dia- 
gram which can be explored with heavy-ion collisions 
at RHIC, SPS, and AGS/SIS300, respectively. In or- 
der to directly compare the QPM with these lattice re- 
sults, we calculate tib from (13) and s from (15, 16) up 
to 0{{fj,B/Tf), where C2,4.6(V) arc obtained from (12), 
co(T) = p{T,iiB = 0)/r4'from (1,2), and the derivatives 
c^(T) are estimated through fine but finite difference ap- 
proximations of the c„(T). 

Besides investigating the impact of different continuum 
extrapolations of Co(T) on the pattern of isentropic tra- 
jectories, we can ask whether the differences observed be- 
tween the parametrizations of cq (T) and C2 (T) can be ab- 
sorbed in such an extrapolation. Note that, even though 
the cutoff dependence of the lattice results is qualitatively 
similar at /Ub = and at ^b 0, no uniform continuum 
extrapolation is expected for the different Taylor expan- 
sion coeflticients [41, 51]. In Fig. 5 we show the raw lattice 
data for co(T) [17] (squares) together with a continuum 
extrapolation (circles) obtained by multiplying the raw 
data for T > Tg by a factor d = 1.1. The correspond- 



ing QPM parametrizations ("fit 1" (dash-dotted) and "fit 
2" (dashed) in the upper panel of Fig. 5) reproduce the 
lattice QCD results impressively well. Nonetheless, the 
corresponding QPM results for C2,4(r) underpredict the 
lattice data, as depicted in the bottom panel of Fig. 5. 
In particular, the pronounced structure in C4(r) at Tc is 
not well reproduced by the QPM fit. If we instead use 
a QPM parametrization that optimally reproduces C2 (T) 
(solid line in the bottom panel of Fig. 5), the correspond- 
ing QPM result for cq{T) ("fit 3" in the upper panel of 
Fig. 5) agrees fairly well with an assumed continuum ex- 
trapolation of the raw lattice data by a factor d = 1.25 
for T > Tc (triangles). 

In Fig. 6, the QPM results for s/ub = 300 and 45 
employing different fits are exhibited together with the 
results of [39]. In the top panel of Fig. 6 we see that 
the lattice results can be fairly well reproduced when 
using simultaneously two separately optimized QPM 
parametrizations for ca(T) and C2(T) (cf. Fig. 1 and 2). 
This approach, however, would give up thermodynamic 
consistency of the model. When using a single consistent 
parametrization for both cq and C2, specifically the one 
shown by the solid lines in Fig. 5 corresponding to "fit 3" , 
the QPM produces the isentropes shown in the bottom 
panel of Fig. 6. (The other two fits shown in Fig. 5 yield 
almost the same isentropic expansion trajectories as "fit 
3".) For large s/ub, i- e. for small net baryon densities, 
differences between the QPM results in the top and bot- 
tom panels of Fig. 6 arc small, although the top fit shows 
a weak structure near Tc which disappears in the selfcon- 
sistent fit shown in the bottom panel. With decreasing 
s/ub the differences between the results from the two 
fitting strategies increase. They are mainly caused by 
differences in the slope of cq{T) which affect the shape 
of s{T)/T^ and translate, for a given isentropic trajec- 
tory, into large variations of hb near Tc(/xb=0) = 175 
MeV while causing only small differences of about 6% 
at large T. In particular, the pronounced structures of 
the isentropic trajectory near the estimated phase bor- 
der are completely lost in the sclfconsistent fit procedure. 
This shows that the pattern of the isentropic expansion 
trajectories is quite sensitive to details of the EoS. For 
instance, when employing co(r) data which were extrap- 
olated to the continuum by multiplication with a factor 
d = 1.25 aX T > Tc while leaving C2Afi{T) unchanged, 
one obtains the isentropic expansion trajectories shown 
by open squares in the bottom panel of Fig. 6 which also 
lack any structure near the phase transition. 

Changing the deconfinement transition temperature 
to Tc = 170 MeV results in a shift of the trajectories 
by about 10% in /ig direction near Tc but has negli- 
gible consequences for T > 1.5Tc. At asymptotically 
large T, where CQfi{T) are essentially flat, the relation 
Es. = i8£a(l]^) holds for small /xs, i. e. Hues of constant 
specific entropy are essentially given by lines of constant 
Hb/T, as is the case in a quark-gluon plasma with per- 
turbatively weak interactions. 

Figure 6 also shows the chemical freeze-out points de- 
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FIG. 5: (Color online) Top panel: co(T) = p{T, fiB=0) /T'^ 
as a function of T/Tc for Nf = 2. Raw lattice QCD data 
from [17] (squares) and guesses for the continuum extrapo- 
lated data obtained by multiplying (for T > = 175 MeV) 
by a factor d = 1.1 (circles) and d — 1.25 (triangles) [17, 40] 
are shown together with the corresponding QPM fits (dashed- 
dotted, dashed, and solid curves, respectively). The QPM 
parameters read B{T^) = 0.31 T*, b = 344.4, A = 2.7, and 
Ts = 0.46Tc for the dashed-dotted line ("fit 1"); they are the 
same as in Fig. 1 for the dashed line ("fit 2"); and the same 
as in Fig. 2 (with B{T^) = 0.61 T^") for the solid line ("fit 
3"). Bottom panel: Corresponding QPM results compared 
with lattice results for C2(T) (squares) and C4(T) (circles) as 
a function of T/Tc with the same line code as in the top panel. 
The horizontal lines indicate the Stefan-Boltzmann values. 



duced from hadron multiplicity data for Au+Au colli- 
sions at = 130 A GeV at RHIC (Tchom = 169 ± 6 
MeV and HB.chcm = 38± 4 MeV [52]) and for 158 A GeV 
Pb-fPb collisions at the CERN SPS (Tchcm = 154.6 ±2.7 
MeV and /ii3,chom = 245.9 ± 10.0 MeV [53]). Note that 
the specific entropies at these freeze-out points as de- 
duced from the statistical model [55] are s /ub — 200 for 
RHIC-130 and s/ub = 30 for SPS-158, i.e. only about 
2/3 of the values corresponding to the QPM fit of the 
QCD lattice data. One should remember, though, that 
the phenomenological values are deduced from experi- 
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FIG. 6: (Color online) Isentropic evolutionary paths. Trian- 
gles and circles indicate Nf = 2 lattice QCD data from [39] 
for s/ns = 300 and 45, respectively. Corresponding QPM 
results are depicted in the upper panel for a mixed fit where 
co{T) and C2{T) were fitted independently (cf. Figs. 1 and 2). 
In the lower panel we show results from "fit 3" from Fig. 5, 
with open squares indicating the corresponding continuum- 
extrapolated lattice results where the raw co{T) lattice data 
were multiplied by a constant factor d — 1.25 at T > Tc [17]. 
Full red squares show chemical freeze-out points deduced in 
[52, 53] from hadron multiplicity data, as summarized in [54]. 



mental data using a complete spectrum of hadronic res- 
onances whereas the lattice simulations were performed 
for only Nf = 2 dynamical quark flavors with not quite 
realistic quark masses. 

Figure 7 shows that along isentropic expansion lines 
the EoS is almost independent of the value of s/ub- 
Accordingly, the speed of sound — dp/de (which 
controls the build-up of hydrodynamic flow) is essen- 
tially independent of the specific entropy. Note that 
whether we employ the mixed fit or the thermodynam- 
ically consistent fits 1, 2 and 3 of Fig. 5 does not sig- 
nificantly affect the EoS along the isentropes; for large 
energy densities e > 30GeV/fm"^ the differences in p{e) 
are less than 2%. 



9 




e [GeV/fm^] 



FIG. 7: (Color online) Lattice QCD data [39] of p as a func- 
tion of e for Nf = 2 along isentropes with s / ub = 300 (tri- 
angles) and 45 (circles), compared with the corresponding 
QPM results (solid blue and dashed black lines, respectively) . 
These two thick lines employ the mixed fit shown in the up- 
per panel of Fig. 6 and are indistinguishable for s/ub =300 
and 45. The thin solid lines show corresponding results for 
the self-consistent "fit 3" from Fig. 5. Again the curves for 
different s/ub are indistinuishable, and also the deviations 
from the mixed fit are minor. 



D. A remark on the QCD critical point 

At a critical point (CP) a first order phase transition 
line terminates and the transition becomes second order. 
QCD with Nf — 2 + 1 dynamical quark flavors with phys- 
ical masses is a theory where such a CP is expected at 
finite T and /is [56-58]. Its precise location is still a mat- 
ter of debate [48, 59-61], but [59] claim Te = 162 MeV 
and ^J,B,E — 360 MeV for the critical values. In the fol- 
lowing, we focus on initial baryon densities ns < 0.5 
fm~'^ which, assuming isentropic expansion with con- 
served s/ub = 250, corresponds to a baryon chemical 
potential /^^(r^ 170 MeV) < 60 MeV. This is sufficiently 
far from the conjectured CP that we should be justified in 
assuming that the EoS is adequately parametrized by our 
QPM for describing bulk thermodynamic properties and 
the hydrodynamical evolution of the hot QCD matter. 



III. EQUATION OF STATE 

In this Section we concentrate on the physical case of 
Nf = 2 + 1 dynamical quark flavors and match the QPM 
fit to the lattice QCD data at temperatures above Tc to a 
realistic hadron resonance gas EoS below Tc. In this way 
we construct an EoS that can be applied to all stages of 
the hydrodynamic expansion of the hot matter created 
in relativistic heavy-ion collisions at RHIC and LHC. We 
focus our attention on the region of small net baryon 
density explored at these colliders. 



A. Pressure as a function of energy density 

Our goal is to arrive at an EoS in the form p(e,nB) 
as needed in hydrodynamic applications. We anchor 
our QPM approach above Tc to lattice QCD simula- 
tions for Nf = 2+1 dynamical quark flavors presented 
in [17, 62, 63] where p{T)/T^ and e{T)/T^ were calcu- 
lated using mgo = 0.4T and = T. Unfortunately, 
Taylor series expansions for non-zero ^b analogous to 
the Nf = 2 case are not available for Nf — 2+1. Ef- 
fects of finite ^B were studied in [64] for Nf = 2 + 1 
by the multi-parameter reweighting method and success- 
fully compared with the quasiparticle model in [37] by 
testing the extrapolation via Eq. (8). We here concen- 
trate on results from lattice QCD simulations employing 
improved actions [17] which strongly reduce lattice dis- 
cretization errors at high temperatures. First, we focus 
on the available data at /xs = and assume that the ex- 
tension to non-zero ^b can be accomplished through the 
QPM without any complications, relying on the success- 
ful test of our model at finite baryon density for Nf = 2 
as reported in the preceding section and earlier publica- 
tions. 

In Fig. 8 we compare the QPM results for the pressure 
p{T)/T^ and entropy density s{T)/T^ with Nf = 2+1 
lattice QCD data where s follows simply from e and p 
through s/T^ = (e+p)/T*. The parametrization found 
at p,B = is now used to obtain the required thermody- 
namic observables at non-zero ub from the full QPM via 
Eqs. (1), (6) and the relation e+p—Ts = ^j-btib, exploit- 
ing the Maxwell relation (8). 

In Fig. 9 we compare the QPM equation of state p{e) 
at ns = with the corresponding lattice QCD result 
deduced from data for p and e at = [17] in the en- 
ergy density domain explored by heavy ion collisions at 
RHIC. The used lattice data [17] were already extrapo- 
lated to the continuum in [63]. In [62, 65] Tc = (173 ±8) 
MeV was found for the deconfinement transition tem- 
perature. Recent analyses [66, 67] have pointed out re- 
maining uncertainties in the extraction of Tc which would 
have to be sorted out by simulations on larger lattices. 
Here, we set the physical scale to Tc — 170 MeV (see 
discussion below). In the transition region the energy 
density e(T) varies by 300% within a temperature in- 
terval of AT « 20 MeV while p{T) rises much more 
slowly (see upper panels in Figs. 8 and 9). This indi- 
cates a rapid but smooth crossover for the phase tran- 
sition from hadronic to quark matter. At large energy 
densities e > 30GeV/fm^ the EoS follows roughly the 
ideal gas relation e = 3p. For the sake of comparison, 
a bag model equation of state describing a gas of mass- 
less non-interacting quarks and gluons with bag constant 
B^/^ = 230 MeV is also shown in Fig. 9 (straight dotted 
line in the top panel). 

As an aside, differences inp(e, ub—O) arising from con- 
sidering different numbers Nf of dynamical quark flavors 
are investigated in the bottom panel of Fig. 9. Compar- 
ing the QPM result for Nf = 2+1 with the result for 
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FIG. 8: (Color online) Comparison of the QPM with lat- 
tice QCD results (symbols) for the scaled pressure p/T*' 
(top panel) and the scaled entropy density s/T^ (bottom 
panel) as a function of T /Tc for Nf = 2 + 1 and = 0. 
The lattice QCD data [63] are already continuum extrap- 
olated. The QPM parameters read A = 7.6, Ts = 0.8Tc, 
h = 348.72 and B{Tc) = 0.52T^'' where = 170 MeV. In the 
top panel, the horizontal line indicates the Stefan-Boltzmann 
value psb/T" = co = (32+21iV/)7rVl80, using iV/=2.5 to 
account for the non-zero strange quark mass. 
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FIG. 9: (Color online) Top paneJ; iV/= 2-1-1 QPM equa- 
tion of state of strongly interacting matter for vanishing 
net baryon density (solid line) compared with Nf = 2+1 
continuum extrapolated lattice QCD data [63] (squares) at 
ns=0. The dotted line represents p{e) for a gas of mass- 
less non-interacting quarks and gluons with a bag constant 
B^^^ = 230 MeV. Bottom panel: QPM EoS for Nf =2 (dashed 
line) employing "fit 2" in Figs. 1 and 5, compared with lattice 
data [63] (squares) and QPM results (solid line) for Nf = 2-1-1, 
in logarithmic representation. 



Nf = 2 (see Fig. 1), the latter exceeds the Nf = 2+1 
result in the transition region (by about 12% at e = 1 
GeV/fm'^). For larger energy densities e > 3 GeV/fm"^ 
the EoS is found to be fairly independent of Nf even 
though at fixed T both p{T) and e(T) are significantly 
smaller for Nf = 2 than for Nf = 2+1 (see Figs. 1, 8). 

B. Baryon density effects 

We turn now to the baryon density dependence of the 
EoS. Since for hydrodynamics the relation p(e, rts) mat- 
ters, we consider the ub dependence of the pressure at 
fixed energy density. Figure 10 shows that significant 
baryon density dependence of the pressure at fixed en- 
ergy density arises only for e < 2 GeV/fm'^. At the small- 



est energy densities considered here, the dependence of 
p on Ub cannot be determined over the entire ub region 
shown since the flow equation (8) for G^(T,iib) has no 
unique solution at large /is for temperatures far below 
the estimated transition temperature Tc{fiB) [38]. How- 
ever, in the family of equations of state that we will con- 
struct and employ in the following, this peculiar feature 
for small e will not occur. Larger baryon densities which 
become relevant at AGS and CERN/SPS energies or the 
future CBM project at the FAIR/SIS300 facility deserve 
separate studies. Under RHIC and LHC conditions fi- 
nite baryon density effects on the equation of state can 
be safely neglected at all energy densities for which the 
QPM model can be applied. 
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FIG. 10: (Color online) Baryon number density dependence 
of the EoS p{e, ub) at constant energy density e as indicated. 
The curves end where the solution of the flow equation (8) is 
not longer unique. 



C. Robustness of the QPM EoS p{e,nB ~ 0) 

We now perform a naive chiral extrapolation of the 
QPM EoS by setting mqo = and ruso = 150 MeV in 
the thermodynamic expressions, leaving all other param- 
eters fixed. The resulting EoS is shown in the top panel 
of Fig. 11. In this procedure a possible dependence of 
the QPM parameters in Eqs. (9), (10) and, especially, of 
the integration constant B(Tc) in Eq. (1) on the quark 
mass parameters ruao is completely neglected. Note that 
in the transition region (e ~ 1 GeV/fm'^) the chirally ex- 
trapolated result exceeds the original QPM equation of 
state (which was fitted to lattice data with unphysical 
quark masses) by approximately 10%. For higher energy 
densities e > 2 GeV/fm^ these quark mass effects are 
seen to be negligible. 

For e < 0.45 GeV/inr^, the fat solid line in the top 
panel of Fig. 11 shows a hadron resonance gas model EoS 
with a physical mass spectrum in chemical equilibrium 
[68]. Obviously, it exceeds both the lattice QCD data and 
their QPM parametrization. The chirally extrapolated 
QPM EoS, on the other hand, approaches and interesects 
the hadron resonance gas EoS. 

Considering p/e as a function of e, we find for the 
lattice-fitted QPM EoS a softest point (p/e),nin = 0.075 
at Be = 0.92 GeV/fm"^. For the chirally extrapolated 
QPM EoS, the softest point moves slightly upward to 
(p/e)min = 0.087 at Cc = l.lGeV/fm'^, in good agree- 
ment with the lattice QCD data which show a softest 
point (p/e)inin = 0.080 at Cc = IGeV/fm^. 

The small differences between the lattice-fitted QPM 
equation of state and its chirally extrapolated version 
for Nf = 2+1 can be further analyzed by studying the 
squared speed of sound . In the middle panel of Fig. 1 1 , 
cj. is shown as a function of T/Tc for both versions of 
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FIG. 11: (Color online) Top panel: QPM EoS for Nf = 2+1 
(solid red) and its chiral extrapolation to physical quark 
masses (dashed blue). Squares show LQCD data for Nf = 
2 + 1 quark flavors with unphysical masses [63] . Middle panel: 
Comparison of the squared speed of sound — dp/de as a 
function of T/Tc from the QPM with lattice QCD data [69] 
(diamonds and triangles) deduced from the Nf = 2 data for 
p{e) in [39]. Differences between the QPM fit to the LQCD 
data (solid red) and its extrapolation to physical quark masses 
(dashed blue) for Nf = 2+1 are almost invisible. Bottom 
panel: Same as middle panel, but plotted as a function of en- 
ergy density e. - In all three panels the solid green line shows 
the hadron resonance gas model EoS "aal" from [68]. 



the QPM EoS and compared with lattice QCD results 
[69]. One sees that, as far as is concerned, the ex- 
trapolation of the QPM to physical quark masses has no 
discernible consequences, and both versions of the QPM 
EOS therefore have identical driving power for collective 
hydrodynamic flow. Hydrodynamically it is thus of no 
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consequence that the available lattice QCD data for the 
EoS were obtained with unphysical quark masses. 

The found EoS is also fairly robust against variations 
in the particular choice of the physical scale Tc- In Fig. 12 
we show p{e) when setting = 160, 170, and 180 MeV, 
respectively, thereby covering the "reasonable range" ad- 
vocated in [62, 65]. For small energy densities and, in 
particular, for large e > 5 GeV/fm'^ the EoS is rather 
independent of the choice of the value for Tc. At inter- 
mediate e, p{e) varies at most by ±20% for AT^ — ±10 
MeV. As discussed below (Section HID), we must any- 
how bridge over this intermediate region when interpolat- 
ing between the QPM and hadron resonance EoS, so this 
weak dependence on the physical scale Tc is irrelevant in 
practice. 
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FIG. 12: (Color online) Dependence of the EoS for Nf = 2-1-1 
on the chosen value of the physical scale Tc- Dashed, full 
and dash-dotted curves correspond to Tc — 160, 170 and 180 
MeV, respectively. Lattice data (squares) from [63]. 

Next we examine variations in p{e, nB~0) arising from 
different continuum extrapolations of the lattice QCD 
data. Considering the various "by hand" continuum ex- 
trapolations of p{T)/T^ shown in Fig. 5 for Nf — 2, the 
resulting EoS are plotted in Fig. 13. Again, some weak 
sensitivity is observed only in the transition region which 
will be bridged over in the next subsection by matching 
the QPM EoS to a realistic hadron resonance gas below 
Tc- The problem discussed in section II B, that different 
optimum QPM parameters are found by fitting the model 
to co(T) or C2(r) (see Figs. 1, 2 and 5), does not matter 
here since the differences in the resulting parametriza- 
tions manifest themselves only weakly in the EoS p{e) 
and are completely negligible for e > 5 GeV/fm'^. In the 
transition region near e ~ 1 GeV/fm'^ the resulting un- 
certainties are of order 20% (see Fig. 13), but again the 
interpolation to the hadronic EoS largely eliminates this 
remaining sensitivity. 

We close this subsection by exploring the robustness of 
the EoS p{e) against variations between different existing 
lattice QCD simulations resulting from present technical 
limitations. In doing so we keep in mind the negligi- 




FIG. 13: (Color online) Dependence of the EoS for N; = 2 
on the employed continuum extrapolation as performed in 
Fig. 5. Dash-dotted, dashed and solid curves correspond to 
the QPM parameterizations of the raw lattice QCD data [17] 
and continuum extrapolations of these data by a factor d = 
1.1 and d = 1.25, respectively. 



bly small baryon density effects in the region ub < 0.5 
fm^'^ pointed out above. In the top panel of Fig. 14 we 
show the available lattice QCD results for p{T)/T^ with 
Nf = 2+1 dynamical quark flavors from three different 
groups [63, 70, 71] and compare them with our QPM 
adjusted individually to each of these data sets. The dif- 
ferences between the data sets reflect the use of different 
lattice actions, lattice spacings, bare quark masses etc. 
As shown in the figure, these differences can be absorbed 
by the QPM through slight variations in the fit parame- 
ters. However, when presenting the lattice results in the 
form of an EoS p(e), they all coincide for e > 5 GeV/fm'^ 
(bottom panel of Fig. 14). [The agreement is excellent 
up to e « 30GeV/fm^ while at higher energy densities 
a small difference of about 6% between the equations of 
state from [63] and [71] begins to become visible.] In this 
region the EoS can be parameterized by p = ae + /3 with 
a = 0.310±0.005 and (3 = -(0.56±0.07) GeV/fm^. This 
robustness of the lattice QCD EoS for e > 5 GeV/fm"^ im- 
plies that it can be considered as stable input for hydro- 
dynamic simulations of heavy-ion collisions, and that the 
equation of state is well constrained at high energy den- 
sities. Our effort to substitute the often used bag model 
EoS above Tc by a realistic QPM EoS which incorporates 
the lattice data seems therefore well justified. 



D. Matching lattice QCD to a hadron resonance 
gas equation of state via the QPM 

In this subsection we will now match the lattice QCD 
EoS at high energy densities with a realistic hadron res- 
onance gas model at low energy densities [72, 73]. Since 
available lattice QCD simulations still employ unrealistic 
quark masses while the hadron gas model builds upon 
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FIG. 14: (Color online) Stability of the QPM EoS fitted to 
lattice QCD results for Nf = 2-\-\. Top panel: The scaled 
pressure p(T)/r* at fj,B~0 from different lattice QCD calcula- 
tions (Ref. [63] (squares), Ref. [70] (diamonds and triangles), 
and Ref. [71] (circles)), together with corresponding QPM 
fits (solid, long-dashed and dash-dotted, and short-dashed 
lines, respectively). The fit parameters are optimized sepa- 
rately in each case, keeping, however, B{Tc) = CSIT^" with 
Tc = 170 MeV in all three parametrizations fixed. Bottom 
panel: The EoS p{e,nB~0) corresponding to the data and 
fits shown in the top panel. 



the measured spectrum of hadronic resonances, we will 
use the QPM to parametrize the lattice QCD EoS and 
extrapolate it to physical quark masses. Such quark mass 
effects matter most at the lower end of the temperature 
range covered by the lattice QCD data which is, however, 
also the region where the transition from the QPM to the 
hadron resonance gas model must be implemented. 

In the vicinity of the phase transition, the conditions 
of the lattice QCD evaluations in Refs. [17, 39] corre- 
spond to a pion mass « 770 MeV. This large pion 
mass reduces the pressure at small energy density below 
that of a realistic hadron resonance gas. Smaller quark 
masses are necessary to properly account for the par- 
tial pressure generated by the light pion modes and their 
remnants in the temperature region around Tc- On the 



other hand, the hadron resonance gas model has been 
shown to be consistent with the QCD lattice data below 
Tc if one appropriately modifies its mass spectrum for 
consistency with the employed lattice parameters [63]. 
We will therefore adopt the hadron resonance gas model 
with physical mass spectrum [72, 73] as an appropriate 
approximation of the hadronic phase [74], and use the 
QPM to parametrize the lattice QCD EoS near and above 
Tc. 

For the hadron resonance gas EoS [72, 73] we use the 
implementation developed for the (2-f l)-dimensional hy- 
drodynamic code package AZHYDRO [68] which pro- 
vides this EoS in tabulated form on a grid in the (e, hb) 
plane. Specifically, we use EoS "aal" from the OSCAR 
website [68] up to ei — 0.45 GeV/fm'^. It describes a 
thermalized, but chemically non-equilibrated hadron res- 
onance gas, with hadron abundance yield ratios fixed at 
all temperatures at their chemical equilibrium values at 
T = Tc = 170 MeV, as found empirically [75] in Au-|-Au 
collisions at RHIC. 

As seen in Fig. 11, the pressure p{e) of the hadron res- 
onance gas EoS does not join smoothly to that of the 
QPM EoS at Tc (i.e. at d = 0.45 GeV/fm^), irrespective 
of whether one uses directly the QPM fit to the lattice 
QCD data with unphysical quark masses (solid red line 
in Fig. 11) or extrapolates the QPM to physical quark 
masses (dashed blue line). A thermodynamically consis- 
tent treatment thus requires a Maxwell like construction, 
equating the two pressures at a common temperature Tc 
and baryon chemical potential ijb- We opt here for a 
slightly different approach which has the advantage of 
allowing a systematic exploration of the effects of details 
(e.g., stiffness or velocity of sound) of the EoS near Tc 
on hydrodynamic flow patterns: We interpolate p{e, ns) 
at fixed baryon density ns linearly between the hadron 
resonance gas ("aal") value at e = ei to its value in 
the QPM at a larger value Cm, keeping ei fixed but let- 
ting the "matching point" value Cm vary. In our pro- 
cedure r(em) > T{ei), so T(e) is also interpolated lin- 
early, as is the baryon chemical potential hb{s) at fixed 
nB- (This is a convenient pragmatic procedure to in- 
terpolate the special tabular forms of the EoS between 
ei and Cm employed below. Complete thermodynamic 
consistency would require involved polynomials for tem- 
perature and chemical potential interpolation. We utilize 
the linearized structures since the hydrodynamical evo- 
lution equations do not explicitly refer to T and fiB hi 
the interpolation region; instead, only p(e,nB) matters.) 

This produces a family of equations of state whose 
members are labelled by the matching point energy den- 
sity Cm- We here explore the range LOGeV/fm"^ < Cm < 
4.0GeY/fm'^ (see Fig. 15). Since the chiral extrapola- 
tion of the QPM fit to physical quark masses signifi- 
cantly affects the EoS p{e) only at energy densities below 
1 GeV/fm'^ (see top panel in Fig. 11), it does not matter 
whether we use for this procedure the direct QPM fit to 
the lattice QCD data or its chiral extrapolation. 

Figure 15 shows the result for four selected Cm val- 
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FIG. 15: (Color online) A family of equations of state for 
Nf =2 + 1, combining our QPM at high energy densities with 
a hadron resonance gas model ("res. gas") in the low energy 
density regime through linear interpolation. We show the 
range of energy densities relevant for collisions at RHIC. The 
solid lines show p{e) for QPM(4.0), QPM(2.0), QPM(1.25), 
and QPM(l.O) (from top to bottom), where the numerical 
label indicates the matching point em in GeV/fm^ On the 
given scale, effects of varying ub between and 0.5 fm~'^ are 
not visible. Lattice QCD data (squares) are from Ref. [63]. 
For comparison a bag model EoS ("bag") with a sharp first 
order phase transition is also shown (dashed line). The bot- 
tom panel zooms in onto the transition region, using a linear 
energy density scale. 



ues. Cm = 1.0, 1.25, 2.0, and 4.0GeV/fm3 (from bot- 
tom to top). For Cm = 3.0 GeV/fm'^ one obtains a curve 
p(e) (not shown) that extrapolates the hadron resonance 
gas with constant slope all the the way to the QPM 
curve. The dashed line in Fig. 15 shows a Maxwell 
construction between the hadron resonance gas and a 
bag model equation of state with = ^; this results 
in a strong first order phase transition with latent heat 



Aeiat = l.lGeV/fm3 ("EoSQ" in [1, 68]). 

Our construction differs from the approach explored in 
[2] where the hadron resonance gas is matched to an ideal 
quark-gluon gas with varying values for the latent heat 
Aeiat- For example, varying the latent heat in EoS Q 
from Aeiat = 0.4GeV/fm3 to 0.8 and l.eGeV/fm^, the 
pressure p{eQ,nB =0) at a typical initial energy density 
eo = 30 GeV/fm"^ for central Au-|-Au collisions at RHIC 
decreases by only 1.4% and 4.3%, respectively, with cor- 
respondingly small changes in the entropy density sq. In 
our approach, however, the entropy density sq at cq is 
given by lattice QGD and significantly (~ 15%) smaller. 
We note that our QPM(l.O) is similar to EOS Q in [1, 68], 
except for the larger latent heat of EoS Q. 



bag 
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FIG. 16: (Color online) Squared speed of sound — dp/de as 
a function of energy density e along an isentropic expansion 
trajectory with s/ub ~ 100, for the EoS family QPM(e,„) 
depicted in Fig. 15. Baryon density effects are not visible on 
the given scale as long as < 0.5 fm~'^. 



Figure 16 shows the corresponding squared speed of 
sound, c^, as a function of energy density e. The lin- 
ear interpolation between the hadron resonance gas at 
e < d = 0.45GeV/fm3 and the QPM at e > Cm leads 
to a region of constant sound speed for ei < e < em- 
This constant increases monotonically with the match- 
ing point value Cm. For Cm = 3GeV/fm'^, the hadron 
resonance gas extrapolates smoothly to the QPM, with 
no "soft region" of small sound speed left over at all. In 
this case the typical phase transition signature of a soft- 
ening of the EoS near Tc is minimized, leading to minimal 
phase transition effects on the development of hydrody- 
namic flow. 
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IV. AZIMUTHAL ANISOTROPY AND 
TRANSVERSE MOMENTUM SPECTRA 

Equipped with our QCD based family of equations of 
state, we can now explore the effects of fine structures 
in the EoS near on the evolution of hydrodynamic 
flow, by computing the transverse momentum spectra 
dN/{dy prdpT d(j)) and elliptic flow V2{pt) for a variety 
of hadron species. To emphasize flow effects, we only 
consider directly emitted hadrons and neglect resonance 
decay distortions. 

In non-central heavy- ion collisions, the initial almond 
shaped cross section of the overlap zone perpendicular 
to the beam direction in coordinate space is converted 
into an azimuthally asymmetric momentum distribution 
due to the appearance of a radially non-symmetric flow 
governed by pressure gradients. Assuming no transverse 
flow at a certain "initial time" tq, at which the hydrody- 
namical expansion stage starts, the azimuthal asymme- 
try is determined by the acting pressure. Therefore, the 
azimuthal asymmetry is an ideal probe of the equation 
of state. In addition, the final anisotropy in the momen- 
tum distribution depends on the rescatterings among the 
particles and serves as measure of the degree of local ther- 
malization. 

The asymmetry is quantified by the harmonic coeffi- 
cients of an expansion of the emitted hadrons transverse 
momentum spectra into a Fourier series in the azimuthal 
emission angle aromid the beam axis relative to the 
reaction plane (which is determined by the direction of 
the impact parameter 6): 

, . . , =7^ ——{l + 2v2{p±,y)cos2cP+...). 

P± dp± dy d(p ZTT p_L dp± dy 

(20) 

The second Fourier coefficient V2{p^_,y) = (cos2(/))p^^y 
is called elliptic flow. We here exploit the 2-1-1 dimen- 
sional relativistic hydrodynamic program package with 
Cooper-Prye freeze-out formalism, AZHYDRO, used in 
Rcfs. [1, 4, 5, 8]. It assumes longitudinally boost- 
invariant expansion a la Bjorken. Clearly, this is ap- 
propriate only near midrapidity ~ 0, but sufficient for 
purposes of our qualitative investigation here. 

Different phenomenological equations of state of 
strongly interacting matter were proposed in previous 
studies [1-5, 8-11, 73], exhibiting either a strong first or- 
der phase transition with different values of latent heats 
[1, 2, 4, 8, 11, 73], a smooth but rapid crossover [5], or no 
phase transition at all [73] . These equations of state dif- 
fer significantly in their high-density regions and softest 
points, and in the speed of sound which controls details of 
the developing flow pattern. Investigating the hydrody- 
namic consequences of different equations of state helps 
to establish benchmarks for tracing specific phase tran- 
sition signatures and distinguishing them from other dy- 
namical features (such as so far poorly explored viscous 
effects). 

We emphasize, however, that we do not attempt here a 



systematic comparison with RHIC data. Previous stud- 
ies [1, 2, 5] have already qualitatively established that 
existing data are best described by an EoS with a phase 
transition or rapid crossover of significant strength (i.e. 
featuring a strong increase of the entropy and energy den- 
sity within a narrow temperature interval) that exhibits 
both a soft part near Tc and a hard part not too far above 
Tc. More quantitative statements about a preference of 
one form of the EoS over another require a discussion 
that goes beyond the pure ideal fluid dynamical approach 
discussed here, due to well-known strong viscous effects 
on the evolution of elliptic flow in the late hadron reso- 
nance gas phase [76] . Studying the effects of EoS varia- 
tions within a more complete framework that allows to 
account for non-ideal fluid behaviour in the very early 
and late stages of the flreball expansion is an important 
task for the future. Staying here within the ideal fluid 
approach, we do note, however, that our discussion im- 
proves over that presented in [5] by employing below Tc 
a chemically non- equilibrated hadron resonance gas EoS 
which correctly reproduces the measured hadron yields, 
irrespective of the selected value for the hydrodynamic 
decoupling temperature. 

A. Top RHIC Energy 

We employ P. Kolb's program package version 0.0 
available from the OSCAR archive [68]. While the 
study presented in [1] shows that at RHIC energies 
(\/i ~ 200 AGeV) most of the finally observed momen- 
tum anisotropy develops before the completion of the 
quark-hadron phase transition, the build-up of elliptic 
flow still occurs mostly in the temperature region where 
the lattice QCD data show significant deviations from an 
ideal quark- gluon gas. It is therefore of interest to inves- 
tigate the effects of these deviations, and of variations of 
the exact shape of the EoS in the transition region, on 
the final elliptic flow in some detail, both at RHIC en- 
ergies, where they are expected to matter, and at higher 
LHC energies where most (although not all [77]) of the 
anisotropic flow will develop before the system enters the 
phase transition region, thus reducing its sensitivity to 
the transition region. 

We fix the initial conditions for top RHIC energy ac- 
cording to [1] 

So = 110fm-^ no = 0.4fm'^ To = 0.6fm/c; (21) 

these parameters describing the initial conditions in the 
fireball center for central {b = 0) Au-I-Au collisions are 
required input for the hydro code [68]. From these ini- 
tial conditions for central collisions the initial profiles for 
non-central collisions are calculated using the Glauber 
model [1]. For our EoS these values translate (indepen- 
dently of the QPM version used) into Cq = 29.8 GeV/fm^, 
Po = 9.4GeV/fm3, and Tq = 357 MeV. [Strictly speak- 
ing, since in the QPM the physical scale is set by Tc, 
varying Tc in the range 170 ± 10 MeV would result in a 
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variation of Cq from 25 to 33 GeV/fm^ when keeping sq 
fixed (such as to maintain the same final charged particle 
multiplicity dN^h/dy cx sqTq). We fix = 170 MeV.] 
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FIG. 17: (Color online) Transverse momentum spectrum 
(left) and elliptic flow coefficient (right) for directly emit- 
ted strange baryons. The symbols represent STAR data [79] 
(see text for details). Solid and dashed curves are for EoS 
QPM(4.0) and the bag model EoS, respectively. 

Our calculations assume zero initial transverse velocity, 
u_LO = at r = tq. In the hadron phase, the Rapp-Kolb 
chemical off-equilibrium EoS [11] is used to account for 
frozen-out chemical reactions. The freeze-out criterion 
is ef.o. = 0.075 GeV/fm^, corresponding to a freeze-out 
temperature of about 100 MeV. All hadrons are assumed 
to freeze out at the same energy density. 

The usual approach when analyzing data is to adjust 
the set of initial and final conditions to keep the trans- 
verse momentum spectra of a given set of hadron species 
fixed, and to then study the variation of vi. Here we in- 
stead illustrate the impact of the EoS by using a fixed set 
of initial and freeze-out parameters. We explore Au+Au 
collisions at a fixed impact parameter h = 5.2 fm, ad- 



justed to best reproduce minimum bias data from the 
STAR collaboration. 

In Fig. 17 we show the transverse momentum spec- 
tra and differential elliptic flow for directly emitted A, 
S, and hyperons. These hadron species do not receive 
large resonance decay contributions, so by comparing the 
results for directly emitted particles with the measured 
spectra one can obtain a reasonable feeling for the level 
of quality of the model description. We show only re- 
sults obtained with the two extreme equations of state, 
QPM(4.0) and the bag model EoS (see Fig. 15). The re- 
sults for QPM(l.O) are very similar to those from the bag 
model EoS, although the latter features a larger latent 
heat. The two remaining equations of state (QPM(1.25) 
and QPM(2.G)) interpolate smoothly between the ex- 
treme cases shown in Fig. 17. 

The left panel in Fig. 17 shows that QPM(4.0) gen- 
erates significantly larger radial flow, resulting in flat- 
ter spectra especially for the heavy hadrons shown 
here. This can be understood from Fig. 15 since this 
EoS does not feature a soft region with small speed of 
sound around T^. Flatter pj- spectra generically result in 
smaller Fourier coefficients V2{pt) [1], but the right panel 
in Fig. 17 shows that for pr < 1.5GeV/c, QPM(4.0) ac- 
tually produces larger V2{pt) than the bag model EoS. 
This implies that QPM(4.0) also produces a larger overall 
momentum anisotropy (i.e. p^-integrated elliptic flow) 
than the bag model EoS, again due to the absence of a 
soft region near Tc- Only at large pr > 2GeV/c, where 
the ideal fluid dynamic picture is known to begin to break 
down [6], does QPM(4.0) give smaller elliptic flow than 
the bag model EoS, as naively expected [1] from the flat- 
ter slope of the single particle pT^-distribution. 

The larger V2{pt) at low pr < 1.5GeV/c from 
QPM(4.G) is not favored by the data. In this sense we 
confirm the qualitative conclusion from earlier studies 
[1, 2, 5] that the data are best described by an EoS with 
a soft region near Tc, followed by a rapid increase of the 
speed of sound Cg above Tc. 

B. LHC estimates 

Predictions for Pb+Pb collisions at the LHC involve a 
certain amount of guesswork about the initial conditions 
at the higher collision energy. We here do not embark 
upon a systematic exploration of varying initial condi- 
tions, as proposed e.g. in Rcfs. [80], but simply guess 
conservatively 

so= 330fm"^ ro = 0.6fm/c, (22) 

keeping all other parameters unchanged. This corre- 
sponds to 3 times larger final multiplicities than mea- 
sured at RHIC. Within the QPM these initial parame- 
ters translate into eo = 127GeV/fm3, po = 42GeV/fm3, 
and To = 515 MeV for the peak values in central Pb-|-Pb 
collisions. We again study collisions at impact param- 
eter b = 5.2 fm, using the Glauber model to calculate 
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FIG. 18: (Color online) Transverse momentum spectrum (left panels) and azimuthal anisotropy (right panels) for pions, kaons 
and protons (upper row) and strange baryons (lower row). Initial conditions according to eq. (22). The spectra show only 
directly emitted hadrons. Solid and dashed curves are for EoS QPM(4.0) and the bag model EoS being similar to QPM(l.O), 
respectively. 



the corresponding initial density profiles from the above 
parameters. 

Again we show results only for the two extreme equa- 
tions of state, QPM(4.0) and the bag model EoS. Gener- 
ally, the pt spectra for LHC initial conditions are flatter 
than for RHIC initial conditions, since the higher initial 
temperature and correspondingly longer fireball lifetime 
results in stronger radial flow. Figure 18 shows that again 
QPM(4.0), which lacks a soft region near Tc, generates 
even larger radial flow (i.e. flatter pT spectra) than the 
bag model EoS (whose results are similar to those ob- 
tained with QPM(l.O)). The radial flow effects are par- 
ticularly strong for the heavy hyperons. 

The overall momentum anisotropy (i.e. the pr- 
integrated elliptic flow) does not increase very much be- 
tween RHIC and LHC [1]. Since the LHC spectra are 
flatter, i.e. have more weight at larger pT than the RHIC 
spectra, the elliptic flow at fixed pT must therefore de- 



crease. This is clearly seen when one compares the right 
panels of Figs. 17 and 18. The decrease is particularly 
strong for the hyperons at low px where the LHC trans- 
verse momentum spectra become extremely flat. 



V. SUMMARY 

We have shown that available lattice QCD calculations 
give converging and robust results for the EoS p{e, ns) in 
the region of large energy density. Baryon density effects 
were shown to be negligibly small for ns < 0.5 fm~^, i.e. 
the EoS relevant for heavy ion collisions at top RHIC 
and LHC energies is the same. In the transition region 
(i.e. for temperatures around Tc) different lattice cal- 
culations still exhibit quantitative differences. The lat- 
tice calculations examined here do not yet join smoothly 
at low energy densities (i.e. at T < Tc) to the hadron 
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resonance gas model EoS with physical mass spectrum. 
While our quasiparticle model covers all considered lat- 
tice QCD equations of state and serves as a reliable tool 
to connect thermodynamic quantities in a thermodynam- 
ically consistent way, it is not obvious that a reliable chi- 
ral extrapolation is feasible by simply rcplac;iiig the quark 
mass parameters employed on the lattice by their phys- 
ical values. If we do so we find significant quark mass 
effects only for energy densities below about 1 GeV/£m^, 
i.e. below the hadronization phase transition. 

In the present paper we therefore assumed as a work- 
ing hypothesis the validity of the hadron resonance gas 
model EoS below Tc (i.e. below an energy density of 
ei = 0.45 GeV/fm'^) and interpolated this EoS linearly 
to the robust high energy density branch from the QPM 
fit to the lattice QCD data. In doing so we arrive at a 
family of equations of state whose members QPM(eni) are 
labeled by the matching point energy density Cm where 
we join the QPM EoS. The resulting equations of state 
QPM(eni) are available in the usual tabulated form on 
the OSCAR website [68] . We find that the uncertain in- 
termediate region, which is bridged over by this interpo- 
lation procedure, has a small but non-negligible impact 
on the evolution of radial and elliptic flow in high energy 
heavy-ion collisions, visible in the transverse momentum 
spectra and elliptic flow coefficients of various (directly 



emitted) hadron species. Existing RHIC data seem to 
favor those members of our family of equations of state 
that exhibit a soft region near Tc followed by a rapid 
rise of the speed of sound towards the ideal gas value 
above Tc. We caution, however, that we did not perform 
a systematic study including simultaneous variations of 
the EoS and initial and final conditions, and that event- 
by-event fluctuations [44-47] or viscous effects [78] may 
wash out differences between different sets of equations 
of state. More quantitative conclusions about the EoS 
require systematic investigations which match the ideal 
fluid description to viscous dynamical models for the very 
early and late stages of the fireball expansion; this is left 
for the future. 
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